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ABSTRACT 

In complex systems with many degrees of freedom such as spin glass and biomolecular 
systems, conventional simulations in canonical ensemble suffer from the quasi-ergodicity 
problem. A simulation in generalized ensemble performs a random walk in potential 
energy space and overcomes this difficulty. From only one simulation run, one can ob- 
tain canonical-ensemble averages of physical quantities as functions of temperature by 
the single-histogram and/or multiple-histogram reweighting techniques. In this article we 
review the generalized-ensemble algorithms. Three well-known methods, namely, multi- 
canonical algorithm, simulated tempering, and replica-exchange method, are described 
first. Both Monte Carlo and molecular dynamics versions of the algorithms are given. We 
then present five new generalized-ensemble algorithms which are extensions of the above 
methods. 

1 INTRODUCTION 

Since the pioneering work of Metropolis and coworkers half a century ago, computer 
simulations have been indispensable means of research in many fields of physical sciences. 
In the field of molecular science, for instance, a number of powerful simulation algorithms 
have been developed (for reviews see, e.g., Refs. 

Canonical fixed-temperature simulations of complex systems such as spin glasses and 
biopolymers are greatly hampered by the multiple- minima problem, or the quasi-ergodicity 
problem. Because simulations at low temperatures tend to get trapped in one of huge 
number of local-minimum-energy states, it is very difficult to obtain accurate canonical 
distributions at low temperatures by conventional Monte Carlo (MC) and molecular dy- 
namics (MD) methods. One way to overcome this multiple-minima problem is to perform 
a simulation in a generalized ensemble where each state is weighted by an artificial, non- 
Boltzmann probability weight factor so that a random walk in potential energy space may 
be realized (for reviews see, e.g., Refs. [S]^|H])- The random walk allows the simulation 
to escape from any energy barrier and to sample much wider configurational space than 
by conventional methods. Monitoring the energy in a single simulation run, one can ob- 
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tain not only the global-minimum-energy state but also canonical ensemble averages as 
functions of temperature by the single-histogram P| and/or multiple-histogram pUl ITT] 
reweighting techniques (an extension of the multiple-histogram method is also referred to 
as weighted histogram analysis method (WHAM) ^5). Besides generalized-ensemble al- 
gorithms, which are usually based on local updates, methods based on non-local updates 
such as cluster algorithms and their generalizations have also been widely used jl2j-jl4j. 
In this article, we focus our discussion on generalized-ensemble algorithms. 

One of the most well-known generalized-ensemble methods is perhaps multicanonical 
algorithm (MUCA) [T^ IT^ (for a review see, e.g., Ref. fZI). (The method is also referred 
to as entropic sampling [TH] , adaptive umbrella sampling jTHj of the potential energy [20] , 
random walk algorithm [211 122], and density of states Monte Carlo [221 • MUCA can also 
be considered as a sophisticated, ideal realization of a class of algorithms called umbrella 
sampling [2^. Also closely related methods are transition matrix methods reviewed in 
Refs. [23IH1-) MUCA and its generalizations have been applied to spin systems (see, e.g., 
Refs. [2ni^[S21)- MUCA was also introduced to the molecular simulation field [HB]- Since 
then MUCA and its generalizations have been extensively used in many applications in 
protein and related systems [S3-[Z1]. Molecular dynamics version of MUCA has also been 
developed [13113120] (see also Refs. [73113] for Langevin dynamics version). MUCA has 
been extended so that flat distributions in other parameters instead of potential energy 
may be obtained [23 12H1 1121 E21 IZHj • Moreover, multidimensional (or multicomponent) 
extensions of MUCA can be found in Refs. [13 113 113 [Hj . 

While a simulation in multicanonical ensemble performs a free ID random walk in 
potential energy space, that in simulated tempering (ST) [73 [77j (the method is also 
referred to as the method of expanded ensemble [73) performs a free random walk in 
temperature space (for a review, see, e.g., Ref. [ZH])- This random walk, in turn, induces 
a random walk in potential energy space and allows the simulation to escape from states of 
energy local minima. ST has also been applied to protein folding problem [7?H 1^ |13 ISO] • 

The generalized-ensemble method is powerful, but in the above two methods the prob- 
ability weight factors are not a priori known and have to be determined by iterations of 
short trial simulations. This process can be non-trivial and very tedius for complex 
systems with many degreees of freedom. Therefore, there have been attempts to ac- 
celerate the convergence of the iterative process for MUCA weight factor determination 
[23 Ha EH Ea E3 123 (see also Refs. [HllHl)- 

In the replica- exchange method (REM) [H3^|H7], the difficulty of weight factor deter- 
mination is greatly alleviated. (A closely related method was independently developed in 
Ref. [HHl- Similar methods in which the same equations are used but emphasis is laid on 
optimizations have been developed [H3 1013] • REM is also referred to as multiple Markov 
chain method [HI] and parallel tempering jTS]- Details of literature about REM and re- 
lated algorithms can be found in recent reviews [^3 13 •) In this method, a number of 
non-interacting copies (or replicas) of the original system at different temperatures are 
simulated independently and simultaneously by the conventional MC or MD method. Ev- 
ery few steps, pairs of replicas are exchanged with a specified transition probability. The 
weight factor is just the product of Boltzmann factors, and so it is essentially known. 

REM has already been used in many applications in protein systems [^3 [Ml IH3 !^5j- 
|lU6j . Other molecular simulation fields have also been studied by this method in various 
ensembles |in7j - [TT3] . Moreover, REM was applied to cluster studies in quantum chem- 
istry field |114j . The details of molecular dynamics algorithm have been worked out for 



REM in Ref. [01] (see also Refs. jOHl HlOj )- This led to a wide application of replica- 
exchange molecular dynamics method in the protein folding problem |115j - P^3] . 

However, REM also has a computational difficulty: As the number of degrees of 
freedom of the system increases, the required number of replicas also greatly increases, 
whereas only a single replica is simulated in MUCA or ST. This demands a lot of computer 
power for complex systems. Our solution to this problem is: Use REM for the weight fac- 
tor determinations of MUCA or ST, which is much simpler than previous iterative methods 
of weight determinations, and then perform a long MUCA or ST production run. The first 
example is the replica- exchange multicanonical algorithm (REMUCA) P7| llU2j . In RE- 
MUCA, a short replica-exchange simulation is performed, and the multicanonical weight 
factor is determined by the multiple-histogram reweighting techniques fOlE]- Another 
example of such a combination is the replica- exchange simulated tempering (REST) pS] . 
In REST, a short replica-exchange simulation is performed, and the simulated tempering 
weight factor is determined by the multiple-histogram reweighting techniques [TUl lllj . 

We have introduced two further extensions of REM, which we refer to as multicanonical 
replica- exchange method (MUCAREM) [0311112] (see also Refs. [1231121]) and simulated 
tempering replica- exchange method (STREM) |128j . In MUCAREM, a replica-exchange 
simulation is performed with a small number of replicas each in multicanonical ensemble 
of different energy ranges. In STREM, on the other hand, a replica-exchange simulation is 
performed with a small number of replicas in "simulated tempering" ensemble of different 
temperature ranges. 

Finally, one is naturally led to a multidimensional (or, multivariable) extension of 
REM, which we refer to as multidimensional replica-exhcange method (MREM) [OH] (see 
also Refs. |129UlU8lll3UllT^ ). A special realization of MREM is replica- exchange umbrella 
sampling (REUS) [OE] and it is particularly useful in free energy calculations. 

In this article, we describe the eight generalized-ensemble algorithms mentioned above. 
Namely, we first review three familiar methods: MUCA, ST, and REM. We then present 
the five new algorithms: REMUCA, REST, MUCAREM, STREM, and MREM (and 
REUS). 



2 GENERALIZED-ENSEMBLE ALGORITHMS 

2.1 Multicanonical Algorithm and Simulated Tempering 

Let us consider a system of atoms of mass {k = 1, ■ ■ ■ , A^) with their coordinate 
vectors and momentum vectors denoted by g = {Qi, ■ ■ ■ , Qn} and p = {Pi, ■ " " jPat}; 
respectively. The Hamiltonian H{q,p) of the system is the sum of the kinetic energy 
K{p) and the potential energy E{q): 

H{q,p)=K{p) + E{q) , (1) 

where 

N 2 

In the canonical ensemble at temperature T each state x = {q,p) with the Hamiltonian 
H{q,p) is weighted by the Boltzmann factor: 



WB{x-T)=exp{-PH{q,p)) , 



(3) 



where the inverse temperature (3 is defined hj (3 = 1/k-BT {k-B is the Boltzmann constant). 
The average kinetic energy at temperature T is then given by 
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Because the coordinates q and momenta p are decoupled in Eq. we can suppress 
the kinetic energy part and can write the Boltzmann factor as 

Wb{x- T) = Wb{E- T) = exp(-/?E) . (5) 

The canonical probability distribution of potential energy Pb{E;T) is then given by the 
product of the density of states n{E) and the Boltzmann weight factor Wb{E; T): 

PBiE;T)(xniE)WBiE;T) . (6) 

Since n{E) is a rapidly increasing function and the Boltzmann factor decreases expo- 
nentially, the canonical ensemble yields a bell-shaped distribution which has a maximum 
around the average energy at temperature T. The conventional MC or MD simulations 
at constant temperature are expected to yield Pb{E;T). A MC simulation based on the 
Metropolis algorithm P is performed with the following transition probability from a 
state X of potential energy £" to a state x' of potential energy E': 

w{x x') = min (^1, ^^1^) = min (1, exp - E)]) . (7) 

A MD simulation, on the other hand, is based on the following Newton equation: 

dE ^ 

Pk= = A, (8) 

where fj^ is the force acting on the A;-th atom {k = 1,---,N). This equation actually 
yields the microcanonical ensemble, and we have to add a thermostat such as Nose- Hoover 
algorithm jlHlf \11V2\ and the constraint method llH4j in order to obtain the canonical 
ensemble. However, in practice, it is very difficult to obtain accurate canonical distri- 
butions of complex systems at low temperatures by conventional MC or MD simulation 
methods. This is because simulations at low temperatures tend to get trapped in one or 
a few of local-minimum-energy states. 

In the multicanonical ensemble [TH| lldj. on the other hand, each state is weighted by 
a non-Boltzmann weight factor Wj^^^lE) (which we refer to as the multicanonical weight 
factor) so that a uniform potential energy distribution Pj^^{E) is obtained: 

Pmn{E) oc n{E)W,nn{E) = coustaut . (9) 

The fiat distribution implies that a free random walk in the potential energy space is real- 
ized in this ensemble. This allows the simulation to escape from any local minimum-energy 
states and to sample the configurational space much more widely than the conventional 
canonical MC or MD methods. 

The definition in Eq. Q implies that the multicanonical weight factor is inversely 
proportional to the density of states, and we can write it as follows: 

WUE) = exp [-(3oEUE; Tq)] = , (10) 

n{E} 



where we have chosen an arbitrary reference temperature, Tq = l/ksPo, and the ''multi- 
canonical potential energt/' is defined by 



E^^{E;To) = kBTolnn{E)=ToS{E) . (11) 

Here, S{E) is the entropy in the microcanonical ensemble. Since the density of states of 
the system is usually unknown, the multicanonical weight factor has to be determined 
numerically by iterations of short preliminary runs ^3 . 

A multicanonical Monte Carlo simulation is performed, for instance, with the usual 
Metropolis criterion P : The transition probability of state x with potential energy E to 
state x' with potential energy E' is given by 

w{x^x') = mm (l. ^'^l^j* ) ^'^^{^'^^) =min(l,exp(-/5oAE^,)) , (12) 
where 

A-Emu = E„iu{E'; To) — E^^{E] Tq) . (13) 

The molecular dynamics algorithm in multicanonical ensemble also naturally follows from 
Eq. (fTUI) . in which the regular constant temperature molecular dynamics simulation (with 
T = Tq) is performed by solving the following modified Newton equation instead of Eq. 
©: iSlllE] 

dE^uiE;To) dE^^^E^To) 

= • ^'^^ 

From Eq. pi|l this equation can be rewritten as 

Pk = T{E) ' ^^"^^ 

where the following thermodynamic relation gives the definition of the "effective temper- 
ature" T{E): 



dS{E) 



= , (16) 

E=E, ^ 



OE 
with 

Ea = <E>TiE.) . (17) 

If the exact multicanonical weight factor WjnniE) is known, one can calculate the 
ensemble averages of any physical quantity A at any temperature T (= I/Zcb/^) as follows: 



< A >T= 



^ A{E)PBiE;T) AiE)n{E)exp{-(3E) 

_E_ _ _E 

Y^Pb{E-T) ~ <E)exp{-(3E) 

E E 



where the density of states is given by (see Eq. |TO|l ) 

The summation instead of integration is used in Eq. ()18|1 . because we often discretize the 
potential energy E with step size e (i? = z = 1, 2, ■ ■ ■). Here, the explicit form of the 



physical quantity A should be known as a function of potential energy E. For instance, 
A{E) = E gives the average potential energy < E >t as a function of temperature, and 
A{E) = l3'^{E- < E >tY gives specific heat. 

In general, the multicanonical weight factor Wyna{E), or the density of states n{E), is 
not a priori known, and one needs its estimator for a numerical simulation. This estimator 
is usually obtained from iterations of short trial multicanonical simulations. The details 
of this process are described, for instance, in Refs. [23 EZ]- However, the iterative process 
can be non-trivial and very tedius for complex systems. 

In practice, it is impossible to obtain the ideal multicanonical weight factor with com- 
pletely uniform potential energy distribution. The question is when to stop the iteration 
for the weight factor determination. Our criterion for a satisfactory weight factor is that 
as long as we do get a random walk in potential energy space, the probability distribution 
Pmu(-E) does not have to be completely fiat with a tolerance of, say, an order of magnitude 
deviation. In such a case, we usually perform with this weight factor a multicanonical 
simulation with high statistics (production run) in order to get even better estimate of the 
density of states. Let A^mu(-E') be the histogram of potential energy distribution P^^{E) 
obtained by this production run. The best estimate of the density of states can then be 
given by the single-histogram reweighting techniques jH] as follows (see the proportionality 
relation in Eq. Q): 

By substituting this quantity into Eq. (fTBj) . one can calculate ensemble averages of phys- 
ical quantity A{E) as a function of temperature. Moreover, ensemble averages of any 
physical quantity A (including those that cannot be expressed as functions of potential 
energy) at any temperature T (= l/k^P) can now be obtained as long as one stores the 
"trajectory" of configurations (and A) from the production run. Namely, we have 

no 

E A{x{k))W-^{E{x{k))) exp [-PE{x{k))] 

< A >T= , (21) 

Y.W-l{E{x{k)))eM-m<m 

k=l 

where x{k) is the configuration at the k-i\i MC (or MD) step and no is the total number of 
configurations stored. Note that when A is a function of E, Eq. (PT|) reduces to Eq. p8|) 
where the density of states is given by Eq. (j2n|l . 

Eqs. fll8|l and ()21|1 or any other equations which involve summations of exponential 
functions often encounter with numerical difficulties such as overflows. These can be 
overcome by using, for instance, the following equation |1351 113b'] : For C = A + B (with 
A > and 5 > 0) we have 



InC =ln 



max(/l, B) 



min(A, B) 
max(A, B) 



[22) 



= max(ln/l. In 5) + In {1 -|- exp [min(lnyl. In B) — max(ln A, In 5)]} . 

We now briefly review the original simulated tempering (ST) method [7^ I77j . In this 
method temperature itself becomes a dynamical variable, and both the conflguration and 
the temperature are updated during the simulation with a weight: 



WsT{E;T) = exp{-PE + a{T)) , 



(23) 



where the function a(T) is chosen so that the probabihty distribution of temperature is 
flat: 

PgT(r) = J dE n{E) Ws-riE] T) = J dE n{E) exp {-pE + a(T)) = constant . (24) 

Hence, in simulated tempering the temperature is sampled uniformly. A free random walk 
in temperature space is realized, which in turn induces a random walk in potential energy 
space and allows the simulation to escape from states of energy local minima. 

In the numerical work we discretize the temperature in M different values, Tm {m = 
1, ■ ■ ■ , M). Without loss of generality we can order the temperature so that Ti < T2 < 
■ ■ ■ < Tm- The lowest temperature Ti should be sufficiently low so that the simulation 
can explore the global-minimum-energy region, and the highest temperature Tm should 
be sufficiently high so that no trapping in an energy-local-minimum state occurs. The 
probability weight factor in Eq. is now written as 

Wst{E; Tm) = exp{-f3mE + a^) , (25) 

where = ci(Tm) (m = 1, ■ ■ ■ , M). Note that from Eqs. (j2H) and we have 

exp(— a^) oc J dE n{E) exp{—PmE) . (26) 

The parameters are therefore "dimensionless" Helmholtz free energy at temperature 
Tm (i.e., the inverse temperature j3m multiplied by the Helmholtz free energy). We re- 
mark that the density of states n{E) (and hence, the multicanonical weight factor) and 
the simulated tempering weight factor am are related by a Laplace transform |13]. The 
knowledge of one implies that of the other, although in numerical work the inverse Laplace 
transform of Eq. ()26p is nontrivial. 

Once the parameters determined and the initial configuration and the initial 

temperature Tm are chosen, a simulated tempering simulation is then realized by alter- 
nately performing the following two steps [7F)| 177] : 

1. A canonical MC or MD simulation at the fixed temperature Tm (based on Eq. ((Zj) 
or Eq. (jH))) is carried out for a certain steps. 

2. The temperature Tm is updated to the neighboring values Tm±i with the configura- 
tion fixed. The transition probability of this temperature-updating process is given 
by the Metropolis criterion (see Eq. p5|l ): 

w{Tm Tm±i) = min (1, exp (-A)) , (27) 

where 

A = {Pm±l - Pm) E - {am±l - ttm) ■ (28) 

Note that in Step 2 we exchange only pairs of neighboring temperatures in order to secure 
sufficiently large acceptance ratio of temperature updates. 

As in multicanonical algorithm, the simulated tempering parameters = o,{Tm) 
(m = 1,---,M) are also determined by iterations of short trial simulations (see, e.g., 
Refs. |78 | I79 | E5] for details). This process can be non-trivial and very tedius for complex 
systems. 



After the optimal simulated tempering weight factor is determined, one performs a 
long simulated tempering run once. The canonical expectation value of a physical quantity 
A at temperature (m = 1, ■ ■ ■ , M) can be calculated by the usual arithmetic mean as 
follows: 

rim 

<A>Tm=—Y.^Mk)) , (29) 

where Xm{k) (/c = 1, ■ ■ ■ ,nm) are the configurations obtained at temperature Tm and rim 
is the total number of measurements made at T = T^. The expectation value at any 
intermediate temperature can also be obtained from Eq. (fTSj) . where the density of states 
is given by the multiple-histogram reweighting techniques ^01 as follows. Let Nm{E) 
and rim be respectively the potential-energy histogram and the total number of samples 
obtained at temperature Tm = l/kB^m ("^ = 1, ■ " " ? M). The best estimate of the density 
of states is then given by (TUl E] 

M 

Y: 9m' Nm{E) 

n{E) = , (30) 

J2 9m exp{fm - (3mE) 
m=l 

where we have for each m (= 1, ■ ■ ■ , M) 

exp(-/„) = n{E) exp(-/5^E) . (31) 

E 

Here, Qm = ^ + 2rm, and Tm is the integrated autocorrelation time at temperature Tm- 
For many systems the quantity Qm can safely be set to be a constant in the reweighting 
formulae [TT], and so we usually set gm = ^■ 

Note that Eqs. (plj) and are solved self-consistently by iteration (TUl E] to obtain 
the density of states n{E) and the dimensionless Helmholtz free energy Namely, we 
can set all the fm ('^ = 1, ■ ■ ■ , M) to, e.g., zero initially. We then use Eq. (j3n|l to obtain 
n{E), which is substituted into Eq. (p?T| to obtain next values of fm, and so on. 

Moreover, ensemble averages of any physical quantity A (including those that cannot 
be expressed as functions of potential energy) at any temperature T (= 1//cb/9) can now 
be obtained from the "trajectory" of configurations of the production run. Namely, we 
first obtain fm {m = 1, ■ ■ ■ , M) by solving Eqs. (j30|l and ()31|) self-consistently, and then 
we have |l()2j 

M nm „-l 

YYMxmik))-, ^- exp [-PE{xm{k))] 

"""^ ''"^ J2 9i'^i^ ^^P ~ l3iE{xm{k))\ 
<^>r= l^m -I ' (32) 

E E -I [-mxm{k))] 

m-1 k=l ^-1^^ _ jj^E{Xm{k))] 

where Xm{,k) (fc = 1, ■ ■ ■ ,nm) are the configurations obtained at temperature Tm- 



2.2 Replica- Exchange Method 

The replica- exchange method (REM) [HSl^IEZI was developed as an extension of simulated 
tempering jHHj (thus it is also referred to as parallel tempering [ZHj) (see, e.g., Ref. [Oil 
for a detailed description of the algorithm). The system for REM consists of M non- 
interacting copies (or, replicas) of the original system in the canonical ensemble at M 
different temperatures (m = 1, ■ ■ ■ , M). We arrange the replicas so that there is always 
exactly one replica at each temperature. Then there exists a one-to-one correspondence 
between replicas and temperatures; the label i {i = 1, ■ ■ ■ , M) for replicas is a permutation 
of the label m (m = 1, ■ ■ ■ , M) for temperatures, and vice versa: 

r i = i(m) = f{m) , .gg. 
1 m = m{i) = f^^{i) , 

where f{m) is a permutation function of m and f~^{i) is its inverse. 

Let X = ^Xi^^^\ ■ ■ ■ , = ■ ■ • , ^^^(m)} stand for a "state" in this general- 

ized ensemble. Each "substate" xj^ is specified by the coordinates g'*' and momenta 
of atoms in replica i at temperature T^: 



x^, 



= . (34) 



Because the replicas are non-interacting, the weight factor for the state X in this 
generalized ensemble is given by the product of Boltzmann factors for each replica (or at 
each temperature): 



M ^ ( M 



W^rem(X) = exp (gl^pl^l) = exp - ^ (gWHl^pW-)]) , (35) 



1=1 ) \ m=l 



where i{m) and m(i) are the permutation functions in Eq. ()33|). 

We now consider exchanging a pair of replicas in the generalized ensemble. Suppose 
we exchange replicas i and j which are at temperatures and T„, respectively: 

Here, i, j, m, and n are related by the permutation functions in Eq. (jSSI), and the exchange 
of replicas introduces a new permutation function /': 

The exchange of replicas can be written in more detail as 



(38) 

where the definitions for pW' and p'-'^' will be given below. We remark that this process is 
equivalent to exchanging a pair of temperatures and Tn for the corresponding replicas 
i and j as follows: 

\ / n \ / m 



In the original implementation of the replica- exchange method (REM) [HSj^lEZ], Monte 
Carlo algorithm was used, and only the coordinates q (and the potential energy function 
E{q)) had to be taken into account. In molecular dynamics algorithm, on the other 
hand, we also have to deal with the momenta p. We proposed the following momentum 
assignment in Eq. (jHEl) (and in Eq. ^) [Hlj: 



p 



p 



Tn 
rr 

T 



(40) 



which we believe is the simplest and the most natural. This assignment means that we 
just rescale uniformly the velocities of all the atoms in the replicas by the square root of 
the ratio of the two temperatures so that the temperature condition in Eq. may be 
satisfied. 

In order for this exchange process to converge towards an equilibrium distribution, 
it is sufficient to impose the detailed balance condition on the transition probability 



VI^rem(^) 



w{X X') 



W^rem(X') 



w{X' X) 



(41) 



Z ' ' z 

where Z is the partition function of the entire system. From Eqs. (^, (0), (jH^j) . pUj) . and 
()41|). we have 



w{X X') 
w{X' X) 



exp |-/3„ 



Kip 



Eiq 



,[i] 



K + E 



where 



A 



exp( 

Pm 
-- iPn 



Tn 

[e 
-A) , 



-Pn[K{p^'^i' 
+ Pn\K(p^^^ 



E{qy^)]} , 



Eiqm -Pn Eiq 



E (g^l) - E (gH)) - /5„ {e (g^l) - E (gt^ 
- Pn) (E (g^l) -E(q 



(42) 

(43) 
(44) 



and i, j, m, and n are related by the permutation functions in Eq. (jH^ before the exchange: 

(45) 



[J =f{n). 

This can be satisfied, for instance, by the usual Metropolis criterion 0: 



w{X ^X')=w 



X, 



[j] 



min (1, exp (—A)) 



(46) 



where in the second expression (i.e., w(a;|^|a;|^^)) we explicitly wrote the pair of replicas 
(and temperatures) to be exchanged. Note that this is exactly the same criterion that 
was originally derived for Monte Carlo algorithm [53]-|87j. 

Without loss of generality we can again assume Ti < T2 < ■ ■ ■ < Tm- A simulation of 
the replica- exchange method (REM) [HHI^IHZ! is then realized by alternately performing 
the following two steps: 



1. Each rephca in canonical ensemble of the fixed temperature is simulated simultaneously 
and independently for a certain MC or MD steps. 

2. A pair of replicas at neighboring temperatures, say and x|^^_,„i, are exchanged 



with the probability w (x]^ ^m+i') Eq. 



Note that in Step 2 we exchange only pairs of replicas corresponding to neighboring tem- 
peratures, because the acceptance ratio of the exchange process decreases exponentially 
with the difference of the two /3's (see Eqs. ()44|1 and ()4fjj) ). Note also that whenever a 
replica exchange is accepted in Step 2, the permutation functions in Eq. are updated. 

The REM simulation is particularly suitable for parallel computers. Because one can 
minimize the amount of information exchanged among nodes, it is best to assign each 
replica to each node (exchanging pairs of temperature values among nodes is much faster 
than exchanging coordinates and momenta). This means that we keep track of the per- 
mutation function m{i] t) = f~^{i] t) in Eq. as a function of MC or MD step t during 
the simulation. After parallel canonical MC or MD simulations for a certain steps (Step 
1), M/2 pairs of replicas corresponding to neighboring temperatures are simulateneously 
exchanged (Step 2), and the pairing is alternated between the two possible choices, i.e., 
(Ti,T2), (T3,r4), ■■■ and {T^,^), {T,,T,),--: 

The major advantage of REM over other generalized-ensemble methods such as mul- 
ticanonical algorithm ^3 ^] and simulated tempering [7^ [77j lies in the fact that the 
weight factor is a priori known (see Eq. while in the latter algorithms the deter- 

mination of the weight factors can be very tedius and time-consuming. A random walk 
in "temperature space" is realized for each replica, which in turn induces a random walk 
in potential energy space. This alleviates the problem of getting trapped in states of 
energy local minima. In REM, however, the number of required replicas increases as the 
system size N increases (according to ^/N) jHS]. This demands a lot of computer power 
for complex systems. 

2.3 Replica- Exchange Multicanonical Algorithm and Replica- 
Exchange Simulated Tempering 

The replica- exchange multicanonical algorithm (REMUCA) ^T7 \ ll()2j overcomes both the 
difficulties of MUCA (the multicanonical weight factor determination is non-trivial) and 
REM (a lot of replicas, or computation time, is required). In REMUCA we first perform 
a short REM simulation (with M replicas) to determine the multicanonical weight factor 
and then perform with this weight factor a regular multicanonical simulation with high 
statistics. The first step is accomplished by the multiple-histogram reweighting techniques 
fOlEI- Let Nm{E) and be respectively the potential-energy histogram and the total 
number of samples obtained at temperature (= l/fce/^m) of the REM run. The density 
of states n{E) is then given by solving Eqs. (jHUj) and (jHTj) self-consistently by iteration. 

Once the estimate of the density of states is obtained, the multicanonical weight 
factor can be directly determined from Eq. (fTTHl (see also Eq. PH) ). Actually, the density 
of states n{E) and the multicanonical potential energy, £'inu(£'; To), thus determined are 
only reliable in the following range: 



El < E < Em , 



(47) 



(4J 



where 

and Ti and Tm are respectively the lowest and the highest temperatures used in the REM 
run. Outside this range we extrapolate the multicanonical potential energy linearly: jHIj 
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E=Ei 



dEjj^uiE; To) 



dE 



(E — Em) + E^^{Em; To 



for E < El, 

for Ei<E< Em, 

for E > Em- 



(49) 



E=Em 



The multicanonical MC and MD runs are then performed respectively with the Metropolis 
criterion of Eq. p2|) and with the modified Newton equation in Eq. in which £^^(E) 
in Eq. is substituted into i?niu(-E'; To). We expect to obtain a fiat potential energy 
distribution in the range of Eq. (jUj). Finally, the results are analyzed by the single- 
histogram reweighting techniques as described in Eq. ^K)^ (and Eq. (fTH|l ). 

Some remarks are now in order. From Eqs. (|TT|). ()16p . (fT7|) . and (jlSj) . Eq. becomes 



' T T 

Y{E - El) + TqS{Ei) = yE + constant , for E < Ei =< E >t^, 

ToS{E) , ' foTEi<E< Em, 

^{E- Em) + ToS{Em) = ^E + constant , for E > Em =< E >t,,. 

^ -l-M -LAI 

(50) 

The Newton equation in Eq. ()14j) is then written as (see Eqs. p5|l . (jlfij) . and ()17|l ) 

(To 
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rp J k ■> 

^ -l-M 



for E < El, 
for Ei<E< Em, 
for E > Em- 



(51) 



Because only the product of inverse temperature f3 and potential energy E enters in the 
Boltzmann factor (see Eq. ©), a rescaling of the potential energy (or force) by a constant, 
say a, can be considered as the rescaling of the temperature by l/a |lHllllUj . Hence, our 
choice of S{°^{E) in Eq. ^ resuhs in a canonical simulation at T = Ti for E < Ei, a 
multicanonical simulation for Ei < E < Em, and a canonical simulation at T = Tm for 
E > Em- Note also that the above arguments are independent of the value of Tq, and we 
will get the same results, regardless of its value. 

For Monte Carlo method, the above statement follows directly from the following 
equation. Namely, our choice of the multicanonical potential energy in Eq. (jl^ gives (by 
substituting Eq. ^ into Eq. (fm|l ) 



iy^,(E) =exp[-/3o4u^(^) 



exp 
1 



-PiE + constant) 



n{E) ' 
exp {-(3mE 



constant) 



for E < El, 

for Ei<E< Em, 

for E > Em- 



(52) 



We now present another effective method of the multicanonical weight factor [2j , which 
is closely related to REMUCA. We first perform a short REM simulation as in REMUCA 
and calculate < E >t as a function of T by the multiple-histogram reweighting tech- 
niques (see Eqs. pO|l and (j3T| ). Let us recall the Newton equation of Eq. (|T5|l and the 
thermodynamic relation of Eqs. and ()17j) . The effective temperature T{E), or the 
derivative ^^^lui^.To) ^ numerically obtained as the inverse function of Eq. ()17p. 

where the average < E >t{e) has been obtained from the results of the REM simulation 
by the multiple-histogram reweighting techniques. Given its derivative, the multicanon- 
ical potential energy can then be obtained by numerical integration (see Eqs. and 



We remark that the same equation was used to obtain the multicanonical weight factor in 
Ref. [S2l! where < E >'r was estimated by simulated annealing instead of REM. Essen- 
tially the same formulation was also recently used in Ref. [Z2] to obtain the multicanonical 
potential energy, where < E >t was calculated by conventional canonical simulations. 

We finally present the new method which we refer to as the replica- exchange simulated 
tempering (REST) [HHl- In this method, just as in REMUCA, we first perform a short 
REM simulation (with M replicas) to determine the simulated tempering weight factor 
and then perform with this weight factor a regular ST simulation with high statistics. 
The first step is accomplished by the multiple-histogram reweighting techniques [TIH ITT] , 
which give the dimensionless Helmholtz free energy fm (see Eqs. flHOj) and (p?T|l ). 

Once the estimate of the dimensionless Helmholtz free energy are obtained, the 
simulated tempering weight factor can be directly determined by using Eq. ()25|) where we 
set ttm = fm (compare Eq. (j26|) with Eq. (j3T| ). A long simulated tempering run is then 
performed with this weight factor. Let Nm{E) and be respectively the potential-energy 
histogram and the total number of samples obtained at temperature (= l/k-BPm) from 
this simulated tempering run. The multiple-histogram reweighting techniques of Eqs. ()30|) 
and pip can be used again to obtain the best estimate of the density of states n{E). The 
expectation value of a physical quantity A at any temperature T (= l/ksP) is then 
calculated from Eq. (fT^ . 

The formulations of REMUCA and REST are simple and straightforward, but the 
numerical improvement is great, because the weight factor determination for MUCA and 
ST becomes very difficult by the usual iterative processes for complex systems. 

2.4 Multicanonical Replica-Exchange Method and Simulated Tem- 
pering Replica-Exchange Method 

In the previous subsection we presented REMUCA, which uses a short REM run for the 
determination of the multicanonical weight factor. Here, we present two modifications 
of REM and refer the new methods as multicanonical replica-exchange method (MU- 
CAREM) p71 llU2j and simulated tempering replica-exchange method (STREM) |128j . 
In MUCAREM the production run is a REM simulation with a few replicas not in the 
canonical ensemble but in the multicanonical ensemble, i.e., different replicas perform 
MUCA simulations with different energy ranges. Likewise in STREM the production run 
is a REM simulation with a few replicas that performs ST simulations with different tem- 



my- m 




(53) 



perature ranges. While MUCA and ST simulations are usually based on local updates, 
a replica-exchange process can be considered to be a global update, and global updates 
enhance the sampling further. 

We first describe MUCAREM. Let A4 be the number of replicas. Here, each replica is 
in one-to-one correspondence not with temperature but with multicanonical weight factors 
of different energy range. Note that because multicanonical simulations cover much wider 
energy ranges than regular canonical simulations, the number of required replicas for the 
production run of MUCAREM is much less than that for the regular REM {Ai -C M). 
The weight factor for this generalized ensemble is now given by (see Eq. i\'d5\\ ) 
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W^mucarem(X) = n ^iu ^^^^ {e {x. 
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where we prepare the multicanonical weight factor (and the density of states) separately 
for m regions (see Eq. (|TIHl ): 



iyi-> (e (xW)) = exp [-/3^4u^ (e (4 
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Here, we have introduced A4 arbitrary reference temperatures = l/ZcB/^m ("^ = 
1,---,A4), but the final results will be independent of the values of T^, as one can 
see from the second equality in Eq. (these arbitrary temperatures are necessary only 
for MD simulations). 

Each multicanonical weight factor W^^{E), or the density of states n^"^^{E), is 
defined as follows. For each m (m = 1,---,A4), we assign a pair of temperatures 



{m} rp{m}^ 
H 



Here, we assume that Tl^"'^ < T^"' and arrange the temperatures so 



■n{m} 



that the neighboring regions covered by the pairs have sufficient overlaps. Without loss 
of generality we can assume T^^^ < ■ ■ ■ < Tjp^-^ and < ■ ■ ■ < . We define the 
following quantities: 
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Suppose that the multicanonical weight factor Wran{E) (or equivalently, the multi- 
canonical potential energy E^^j^{E; Tq) in Eq. (fTT|) ) has been obtained as in REMUCA or 
by any other methods in the entire energy range of interest (-^l^^ < E < Ej;^^). We then 
have for each m (m = 1, ■ ■ ■ , A^) the following multicanonical potential energies (see Eq. 

(gni)): inn 
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(57) 

Finally, a MUCAREM simulation is realized by alternately performing the following 
two steps. 



1. Each replica of the fixed multicanonical ensemble is simulated simultaneously and 
independently for a certain MC or MD steps. 

2. A pair of replicas, say i and j, which are in neighboring multicanonical ensembles, 
say m-th and (m+l)-th, respectively, are exchanged: X = |- ■ ■ , a:^, ■ ■ ■ , x\;ll_^_i, ■ ■ - j — > 

X' = {■ ■ ■ , x^m^ ■ ■ ■ ' ^m+i5 ■ ■ ■}• The transition probabihty of this replica exchange 
is given by the Metropolis criterion: 

w(X ^ X') = min(l,exp(-A)) , (58) 

where we now have (see Eq. fl4Hj) ) [HZj 

A = Pm {Sj-^ {E - {E {^^^^ {e (gt^')) " (^'))} " 

(59) 

Here E (^gt*^ j and E (^q^^^^ are the potential energy of the i-th replica and the j-th 
replica, respectively. 

Note that in Eq. ()59j) we need to newly evaluate the multicanonical potential energy, 
Sj^^E{q^^^)) and Sj^+^^Eiq^"^)), because Sj^^E) and S^JiE) are, in general, different 
functions for m 7^ n. 

In this algorithm, the m-th multicanonical ensemble actually results in a canonical 
simulation at T = T^"^^ for E < Ej\^\ a multicanonical simulation for Ej\^^^ ^ E < E^\ 
and a canonical simulation at T = for E > Ejj^\ while the replica-exchange process 
samples states of the whole energy range {Ei'^ <E< 4-^^ 

For obtaining the canonical distributions at any intermediate temperature T, the 
multiple-histogram reweighting techniques [TUl E] are again used. Let Nm{E) and 
be respectively the potential-energy histogram and the total number of samples obtained 
with the multicanonical weight factor W^^{E) (m = 1, ■ ■ ■, A1). The expectation value 
of a physical quantity A at any temperature T (= I/Zcb/?) is then obtained from Eq. (fTH|) . 
where the best estimate of the density of states is obtained by solving the WHAM equa- 
tions, which now read jHZI 

M M 

Y: g-'N^iE) NUE) 

n{E) = ^ = ^ ^-^ , (60) 

Y g-'nmeMfm)Wi^\E) Y a^n^ exp (/„ - /3^^4u^(^)) 

m=l m=l 

and for each m (= 1, ■ ■ ■ , A^) 

exp(-/„) = Y.<E)Wl:^\E) = Y: n{E) exp {-P^SL^^He)) . (61) 

E E 

Note that W^^^{E) is used instead of the Boltzmann factor exp(— /Sm-E) in Eqs. (jHUj) and 

Moreover, ensemble averages of any physical quantity A (including those that cannot 
be expressed as functions of potential energy) at any temperature T (= I/Zcb/^) can now 
be obtained from the "trajectory" of configurations of the production run. Namely, we 



first obtain (m = 1, ■ ■ ■ ,^A) by solving Eqs. (|Hn|) and (jHT|l self-consistently, and then 
we have |lU2j 

M rim „-l 

y: e AixUk))^ ^- exp [-m^mim 

< A >T= n — ^ , (62) 

M rim „-l ' ^ ^ 

E E M i-m^mim 

i=l 

where the trajectories Xm{k) {k = 1, - ■ ■ , rim) are taken from each multicanonical simula- 
tion with the multicanonical weight factor W^^{E) (m = 1, ■ ■ ■ ,Ai) separately. 

As seen above, both REMUCA and MUCAREM can be used to obtain the multi- 
canonical weight factor, or the density of states, for the entire potential energy range of 
interest. For complex systems, however, a single REMUCA or MUCAREM simulation is 
often insufficient. In such cases we can iterate MUCA (in REMUCA) and/or MUCAREM 
simulations in which the estimate of the multicanonical weight factor is updated by the 
single- and/or multiple-histogram reweighting techniques, respectively. 

To be more specific, this iterative process can be summarized as follows. The RE- 
MUCA production run corresponds to a MUCA simulation with the weight factor Wmu{E). 
The new estimate of the density of states can be obtained by the single-histogram reweight- 
ing techniques of Eq. (j20|l . On the other hand, from the MUCAREM production run, 
the improved density of states can be obtained by the multiple-histogram reweighting 
techniques of Eqs. (jHUj) and (jHH). 

The improved density of states thus obtained leads to a new multicanonical weight 
factor (see Eq. (fTUI) ). The next iteration can be either a MUCA production run (as in 
REMUCA) or MUCAREM production run. The results of this production run may yield 
an optimal multicanonical weight factor that yields a sufficiently flat energy distribution 
for the entire energy range of interest. If not, we can repeat the above process by obtaining 
the third estimate of the multicanonical weight factor either by a MUCA production run 
(as in REMUCA) or by a MUCAREM production run, and so on. 

We remark that as the estimate of the multicanonical weight factor becomes more 
accurate, one is required to have a less number of replicas for a successful MUCAREM 
simulation, because each replica will have a flat energy distribution for a wider energy 
range. Hence, for a large, complex system, it is often more efficient to first try MU- 
CAREM and iteratively reduce the number of replicas so that eventually one needs only 
one or a few replicas (instead of trying REMUCA directly from the beginning and iterat- 
ing MUCA simulations). 

We now describe the simulated tempering replica-exchange method (STREM) |128j . 
Suppose that the simulated tempering weight factor Wst{E; Tn) (or equivalently, the di- 
mensionless Helmholtz free energy a„ in Eq. ()25p) has been obtained as in REST or by 
any other methods in the entire temperature range of interest (Ti < T„ < Tm)- We devide 
the overlapping temperature ranges into regions {M. <^ M). Suppose each tempera- 
ture range m has Mm temperatures: T^^^ (/c = 1, ■ ■ ■ , A/'m) for m = 1, ■ ■ ■ , M.. We assign 



each temperature range to a replica; each rephca i is in one-to-one correspondence with 
a different temperature range m of ST run, where Ti"^^ < T^!™^ < Tjj^ (/c = 1, ■ ■ ■ ,J\fm)- 
We then introduce the rephca-exchange process between neighboring temperature ranges. 
This works when we allow sufficient overlaps between the temperature regions. 

A STREM simulation is then realized by alternately performing the following two 
steps. [I2H1 

1. Each replica performs a ST simulation within the fixed temperature range simultaneously 
and independently for a certain MC or MD steps. 

2. A pair of replicas, say i and j, which are at, say T = T^™^ and T = T^"^'^^\ in neigh- 
boring temperature ranges, say m-th and (m + l)-th, respectively, are exchanged: 
X = ^- ■ ■ , x^l\ ■ ■ ■ , x^f, ■ ■ - j — y X' = {■ ■ ■ , x^k, ■ ■ ■ 5 ■ ■}• The transition prob- 
ability of this replica exchange is given by the Metropolis criterion: 

w(X X') = min(l,exp(-A)) , (63) 

where 

While in MUCAREM each replica performs a random walk in multicanonical ensemble 
of finite energy range, in STREM each replica performs a random walk by simulated 
tempering of finite temperature range. These "local" random walks are made "global" to 
cover the entire energy range of interest by the replica-exchange process. 

2.5 Multidimensional Replica-Exchange Method 

We now present our multidimensional extension of REM, which we refer to as multidi- 
mensional replica- exchange method (MREM) j^ni- The crucial observation that led to the 
new algorithm is: As long as we have M non-interacting replicas of the original system, 
the Hamiltonian H{q,p) of the system does not have to be identical among the replicas 
and it can depend on a parameter with different parameter values for different replicas. 
Namely, we can write the Hamiltonian for the i-th replica at temperature as 

if,„(gW,/]) = A'(/]) + EA„(gW) , (65) 

where the potential energy Ex^ depends on a parameter and can be written as 

E,Jq^^)=Eo{q^^) + X^V{q^^) . (66) 

This expression for the potential energy is often used in simulations. For instance, in 
umbrella samphng [211, ^o{q) and V{q) can be respectively taken as the original potential 
energy and the "biasing" potential energy with the coupling parameter A^. In simulations 
of spin systems, on the other hand, Eo^q) and V{q) (here, q stands for spins) can be 
respectively considered as the zero-field term and the magnetization term coupled with 
the external field A^- 

While replica i and temperature are in one-to-one correspondence in the original 
REM, replica i and "parameter set" A^ = (T^, A^) are in one-to-one correspondence 
in the new algorithm. Hence, the present algorithm can be considered as a multidimen- 
sional extension of the original replica-exchange method where the "parameter space" is 



one-dimensional (i.e., = T^). Because the replicas are non-interacting, the weight 
factor for the state X in this new generalized ensemble is again given by the product of 
Boltzmann factors for each replica (see Eq. (jH^j) ): 

r M 

V^^mrem(^) = exp <^ -Y,Pm{i)Hmii) {q^^,P^^) 

= exp |- 5] P^H^ ^qHn^)]^pl^im)]^^ I ^ 

I 771=1 J 

where i{m) and m{i) are the permutation functions in Eq. Then the same derivation 

that led to the original replica-exchange criterion follows, and the transition probability 
of replica exchange is given by Eq. PUI) . where we now have (see Eq. fj43p ) [HE] 

Here, and are the total potential energies (see Eq. (|66|l ). Note that we need to 
newly evaluate the potential energy for exchanged coordinates, E\^{q^^) and -E'a„(5'^'0, 
because and E\^ are in general different functions. 

For obtaining the canonical distributions, the multiple-histogram reweighting tech- 
niques [ini E] are particularly suitable. Suppose we have made a single run of the 
present replica-exchange simulation with M replicas that correspond to M different pa- 
rameter sets Km = (Tm, Am) ("2 = 1, ■ ■ ■ , M). Let Nm{Eo, V) and Um be respectively the 
potential-energy histogram and the total number of samples obtained for the m-th pa- 
rameter set Am- The WHAM equations that yield the canonical probability distribution 
Pt^\{Eq,V) = n{EQ,V) exp{—PEx) with any potential-energy parameter value A at any 
temperature T = l/k-^P are then given by [211 

M 

g-' NmiEo,V) 
n{Eo, V) = 

Y 9m exp {fm - f^mExJ 

m=l 

and for each m (= 1, ■ ■ ■ , M) 

exp(-/m) = E ^(Eo, V) exp {-f^mExJ . (70) 

Eo,V 

Here, n{EQ,V) is the generalized density of states. Note that n{EQ,V) is independent 
of the parameter sets A^ = (T^, A^) {m = ^,---,M). The density of states n{EQ,V) 
and the "dimensionless" Helmholtz free energy fm in Eqs. fj69| and (|Tn|l are solved self- 
consistently by iteration. 

Incidentally, these formulations of MREM give multidimensional extensions of RE- 
MUCA p7l llU2j and REST ^Hl. In the former, we obtain uniform distributions both 
in Eq and V, whereas in the latter, the parameter sets A^ become dynamical variables 
and a uniform distribution in those parameters will be obtained. Namely, after a short 
MREM simulation, we can use the multiple-histogram reweighting techniques of Eqs. 
(jnni) and ()7()|) to obtain n{EQ, V) and fm- Hence, we can determine the multidimensional 
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multicanonical weight factor Wjnu{Eo,V) and the multidimensional simulated tempering 
weight factor Wst{Eo, V\ A^)- The former is given by 

WUE^,V)= ^ , (71) 
n(Eo, V) 

and the latter is given by (see Eq. (f^ ) 

W^st(^o, y\ = exp i-/3^E^^ + U) • (72) 

We can use MREM for free energy calculations. We first describe the free-energy 
perturbation case. The potential energy is given by 

Exiq) = EM) + A {Ep{q) - Ei{q)) , (73) 

where Ej and Ep are the potential energy for a "wild-type" molecule and a "mutated" 
molecule, respectively. Note that this equation has the same form as Eq. ()fif)|l . 

Our replica-exchange simulation is performed for M replicas with M different values 
of the parameters = (Tm,Xm)- Since Ex=o{q) = Ej{q) and Ex=i{q) = Ep{q), we 
should choose enough values distributed in the range between and 1 so that we 
may have sufficient acceptance of replica exchange. From the simulation, M histograms 
Nm{Ej, Ep — Ei), or equivalently Nm{Ej, Ep), are obtained. The Helmholtz free energy 
difference of "mutation" at temperature T (= l/kB/3), AF = Fx=i — -Fa=o, can then be 
calculated from 

Pt,x=i{Ei,Ep) 

exp(-/?AF) = ^ = ^ . , (74) 

^T,A=0 2^ }^T,\=o{El,Ep) 

Ej,Ep 

where Pt,\{Ei, Ep) = n{Ei, Ep) exp {—(3E\) are obtained from the WHAM equations of 
Eqs. dnni) and dZni). 

We now describe another free energy calculations based on MREM applied to umbrella 
sampling [21], which we refer to as replica- exchange umbrella sampling (REUS). The 
potential energy is a generalization of Eq. and is given by 

E^{q)=Eo{q)+j:X^'^Ve{q), (75) 

e=i 

where Eo{q) is the original unbiased potential, Ve{q) {i = 1, ■ ■ ■ , L) are the biasing (um- 
brella) potentials, and A^^'' are the corresponding coupling constants (A = (A^^-*, ■ ■ ■ , A*-^^)). 
Introducing a "reaction coordinate" ^, the umbrella potentials are usually written as har- 
monic restraints: 

Veiq) = h{aq)-di)' , (£ = 1,---,L) , (76) 

where dg are the midpoints and kg are the strengths of the restraining potentials. We 
prepare M replicas with M different values of the parameters = (T^, Am), and the 
replica-exchange simulation is performed. Since the umbrella potentials V^(g) in Eq. (f7H|) 
are all functions of the reaction coordinate C, only, we can take the histogram Nm{EQ,C,) 



instead of Nm{Eo, Vi,--- , 14). The WHAM equations of Eqs. ^ and ^ can then be 
written as jnHl 

M 

niEo, = — (77) 

J2 9m exp - PmE)^^ 

m=l 

and for each m (= 1, ■ ■ ■ , M) 

exp(-/„) = J2 exp {-(3mEx^) • (78) 

The expectation value of a physical quantity A with any potential-energy parameter value 
A at any temperature T (= l/fce/?) is now given by 

Y^A{E,,i)P^yS^E,,i) 

<^>t\ = ^^'^ ; ^ > (79) 

' E^tA(^o,0 



where ^{E^, ^) = n{EQ,^) exp \^—(3E^j is obtained from the WHAM equations of 
Eqs. dZZj) and (ff8|l . 

The potential of mean force (PMF), or free energy as a function of the reaction coor- 
dinate, of the original, unbiased system at temperature T is given by 



E^T,A={0}(^0,0 

Eo 



(80) 



where {0} = (0, ■■■,0). 

We now present two examples of realization of REUS. In the first example, we use only 
one temperature, T, and L umbrella potentials. We prepare replicas so that the potential 
energy for each replica includes exactly one umbrella potential (here, we have M = L). 
Namely, in Eq. (f75j) for A = A^ we set 

^li^ = k,m , (81) 

where 5k^i is Kronecker's delta function, and we have 

i?;^^(gH) = Eo(g[^l) + V;n,(g[*l) . (82) 

We exchange replicas corresponding to "neighboring" umbrella potentials, Vm and Vm+i- 
The acceptance criterion for replica exchange is given by Eq. ()46p. where Eq. ()68|) now 
reads (with the fixed inverse temperature /? = 1 / k-eT) 



A = /3 (Kn (g'^'l) - Vm. (gf^l) - (g'^'^) + K^+i (g^^^)) , (83) 

where replica i and j respectively have umbrella potentials Vm and Kn+i before the ex- 
change. 



In the second example, we prepare Nt temperatures and L umbrella potentials, which 
makes the total number of replicas M = x L. We can introduce the following re- 
labeling for the parameters that characterize the replicas: 

Am = (Tm, Am) > Ajj = (T/, Aj) . , . 

(m = l,---,M) {I = 1,---,Nt, J=1,---,L) ' 

The potential energy is given by Eq. ()82p with the replacement: m ^ J. We perform the 
following replica-exchange processes alternately: 

1. Exchange pairs of replicas corresponding to neighboring temperatures, T/ and T/+i 
(i.e., exchange replicas i and j that respectively correspond to parameters A/^j and 
A/_|_i j). (We refer to this process as T-exchange.) 

2. Exchange pairs of replicas corresponding to "neighboring" umbrella potentials, Vj 
and Vj+i (i.e., exchange replicas i and j that respectively correspond to parameters 
A/ J and A/ j+i). (We refer to this process as A-exchange.) 

The acceptance criterion for these replica exchanges is given by Eq. ()4fij) . where Eq. ()68|) 
now reads [02] 

A = {(3j- {Eo + Vj (gt^l) - Eo (g^) - Vj (g^)) , (85) 

for T-exchange, and 

A = {Vj - Vj (gW) - Vj^i (gt^l) + Vj+, (g^)) , (86) 

for A-exchange. By this procedure, the random walk in the reaction coordinate space as 
well as in the temperature space can be realized. 



3 CONCLUSIONS 

In this article we have reviewed uses of generalized-ensemble algorithms for both Monte 
Carlo simulations and molecular dynamics simulations. A simulation in generalized en- 
semble realizes a random walk in potential energy space, alleviating the multiple-minima 
problem that is a common difficulty in simulations of complex systems with many degrees 
of freedom. 

Detailed formulations of the three well-known generalized-ensemble algorithms, namely, 
multicaonical algorithm (MUCA), simulated tempering (ST), and replica-exchange method 
(REM), were given. 

We then introduced five new generalized-ensemble algorithms that combine the merits 
of the above three methods. We refer to these methods as replica-exchange multicanoni- 
cal algorithm (REMUCA), replica-exchange simulated tempering (REST), multicanonical 
replica-exchange method (MUCAREM), simulated tempering replica-exchange method 
(STREM), and multidimensional replica-exchange method (MREM), the last of which 
also led to replica-exchange umbrella sampling (REUS). 

The question is then which method is the most recommended. We have recently 
studied the effectiveness of MUCA, REM, REMUCA, and MUCAREM in the protein 
folding problem |l()2j . Our criterion for the effectiveness was how many times the random 



walk cycles between the high-energy region and low-energy region are realized within a 
fixed number of total MC (or MD) steps. We found that once the optimal MUCA weight 
factor is obtained, MUCA (and REMUCA) is the most effective (i.e., has the most number 
of random walk cycles), and REM is the least jin2j . We also found that once the optimal 
ST weight factor is obtained, ST (and REST) has more random walk cycles than REM 
|98| I128j . Moreover, we compared the efficiency of Berg's recursion |HS] , Wang-Landau 
method EH 1221, and REMUCA/MUCAREM as methods for the multicanonical weight 
factor determination in two-dimensional 10-state Potts model and found that the three 
methods are about equal in efficiency |137j - p!^ . 

Hence, the answer to the above question will depend on how much time one is willing 
to (or forced to) spend in order to determine the MUCA or ST weight factors. Given a 
problem, the first choice is REM because of its simplicity (no weight factor determination 
is required). If REM turns out to be insufficient or too much time-consuming (like the 
case with first-order phase transitions), then other more powerful algorithms such as those 
presented in the present article are recommended. 
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